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^ I Abstract 

\o ■ 

i A class of finite difference schemes for solving a fractional anti-diffusive equation, recently pro- 

posed by Andrew C. Fowler to describe the dynamics of dunes, is considered. Their linear stability 
is analyzed using the standard Von Neumann analysis: stabiUty criteria are found and checked nu- 
merically. Moreover, we investigate the consistency and convergence of these schemes. 
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\Q ■ 1 Introduction 

oo " 

■ Partial Differential Equations with nonlocal or fractional operators are widely used to model scientific 

■ problems in mechanics, physics, signal processing, see for example ||3l and references therein. We 
^ ' consider in this chapter a nonlocal conservation law which appears in the formation and dynamics of 

sand structures such as dunes and ripples |j7] [TT]. Since it is generally impossible to obtain analyti- 
cal solutions of these nonlocal models, one must rely on numerical solutions. In the last few decades, 
significant advances in numerical analysis and computational implementation of numerical methods for 
nonlocal/fractional PDEs have been made. For instance, ||6] propose a finite volume method to approxi- 
mate the solutions of a fractal scalar conservation law, that is to say a conservation law regularized by a 
diffusive fractional power of the Laplacian operator and llT3l [T5l use finite difference methods to approx- 
imate fractional diffusive equations. 

In this chapter, we develop the basic numeiical analysis of the following evolution equation proposed by 
A.C. Fowler (see Q, ||8l and ||9l for more details) to study the nonlinear dune formation: 



dtu{t,x) + [ — ] {t,x) + r]I[u{t,-)]{x) -ed^^u{t,x) = 0, (1) 
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where u = u{t, x) represents the dune height and X is a nonlocal operator defined as follows: for any 
Schwartz function ip G 5(M) and any x G M, 



i[v]{x) := / icr3<^"(x-ode (2) 

Jo 

The second and fourth terms of equation ([T|l correspond to the nonlinear and dissipative terms respec- 
tively, while the third term is the nonlocal term, which is anti-dissipative as we will show later on. The 
positive parameters e (resp. t]) quantify the amount of local diffusion (resp. nonlocal anti-diffusion). 

Remark 1. For causal functions (i.e. (p{x) = for x < 0), this operator is, up to a multiplicative 
constant, the Riemann-Liouville integral which is defined as follows: 

1 "(^_^) ^-2/3 ^4/3 



r(|) io ^ ' dx^/^ 

with r the Euler function. 

Many numerical methods for the evaluation of fractional order integrals and the solution of fractional 
order equations are proposed in the literature. Usually, time and spatial fractional derivatives are consid- 
ered: we refer for instance to ll5] [T2][T5l . 

In our case, the integral operator X can be seen as a fractional power of order 2/3 of the Laplacian with 
the bad sign. Indeed, it has been proved that X has the following Fourier transform lH] : 

T{X[^]m=Mi)^v{i), (4) 

where V'xlO = -ax\i\^i +ibxi\i\^i with ax = 2 7r2r(|), 5x = 2 7r2 V3r(|) and J" denotes the Fourier 
transform defined for / G L^{^) by: for all ^ G M 

-^/(e)= / e-2--«/(x)dx. 
Formula dU stems from the following integral formula lH] : 

2:M(x) = - dz. (5) 

Finally, equation ([T|) involves two antagonistic terms: the anti-diffusive operator X which creates insta- 
bilities and the diffusion operator —d^.^ which controls these perturbations. 

Recently, some theoretical results regarding the Fowler model ([T]l have been obtained, namely, existence 
of travelling-waves, the global well-posedness, the failure of the maximum principle and the instability 
of constant solutions |[T]l2lSl- The last two results are a consequence of the non-positivity of the kernel 
K of X - dl^ defined for t > and x G M by 

K{t,-){x) :=^-i(e-*(^"'l-l^+^i{-)))(^). (6) 

These two "bad properties" show that the discrete problem must be handled with care. Indeed, for mono- 
tone models, a classical way to get numerical stability criteria for explicit scheme is to ensure that the 
approximated problem satisfies the discrete maximum principle, which cannot be true for Equation ([TJ. 
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In m, some numerical results regarding this equation have been obtained using an explicit finite differ- 
ence scheme but the detailed numerical study was not performed. Hence, in this chapter, we would like 
to go one step further investigating the numerical stability, consistency and convergence of a class of 
explicit finite difference schemes approximating the Fowler equation. 

The numerical stability is specially interesting here because the growth of the solution depends on fre- 
quencies and time. Hence, the notion of ^-stability, also called strong stability, is not suitable nor 
desirable. In the literature, some authors use another definition of stability, less restrictive than the A- 
stability: the C-stability. This is an abbreviation for convergence stability and is linked with stability in 
the Lax-Richtmyer sense. In this definition, a numerical scheme is stable for the norm 1 1 • 1 1 if for all 
T > 0, there exists a constant K(T) > independent of the time and space steps 6x, 6t such that for all 
initial data uq 

||n"|| <i^(T)||^zo||, VO<?i<^, 

where represents the approximated solution at the time tn = n6t. This definition allows the solution 
to grow with time, which is the case for example for the equation ut — Uxx = cu. For the L^-stability, a 
simple way to prove the numerical stability and specially to get stability criteria is Fourier analysis, see 
Section [3l Hence, considering the C-stability, the Von Neumann condition is written as 

3C > 0,36t* > 0, such that V(^t e]0,6t*]-yk G Z 

\g{k)\ < 1 + C6t, (7) 

where g is the discrete amplification factor, k the wave number and C is a positive constant independent 
of 5x and St. If C = 0, the Von Neumann condition coincides with the A-stability. 
As we will see later in Section |2] the amplification of solutions of the Fowler equation also depends on 
frequencies: low frequencies are slowly amplified whereas the high frequencies are dampened. Hence, 
the notion of C-stability is not adapted for this model because it considers only the amplification due to 
time. To take into account this phenomenon, the "constant" C introduced in the Von Neumann condition 
(|7]l should also depend on the space step in order to be able to control the amplification w.rt. different 
frequencies and this is not possible for a constant, by definition. Since high frequencies are usually 
responsible of numerical instabilities, we are going to focus our attention on them. Thereafter, the idea 
is to exhibit numerical stabihty conditions to ensure the validity of simulations. We then seek numerical 
stability criteria such that the amplification factor satisfies: 

V|A:| > A:o,|5(A:)| < 1, (8) 

where /cq is some threshold frequency. To ensure this inequality, we will exhibit two sufficient condi- 
tions. The first one is rather unusual: it imposes to the space step 6x to be smaller than a given positive 
constant which depends on the ratio e/r] of local diffusion to non-local anti-diffusion. We will in fact 
check numerically that this condition is not necessary. The second one looks more familiar. It is a classi- 
cal CFL-type condition modified by ai]6t/ 5x^/'^ term, which stems from the nonlocal operator We will 
see in the numerical simulations that this condition is both necessary and sufficient to ensure numerical 
stability. 

For a comprehensive study, we also carry out an enor analysis: we compute the truncation and phase 
errors of several finite difference schemes. We finally investigate the convergence of these schemes. 

The remaining of this chapter is organized as follows: in the next section, we present finite difference 
schemes with some discrete version for the fractional derivative and we study the continuous amplifica- 
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tion factor of the linearized Fowler model. Sections |3]and|4]are, respectively, devoted to the stability and 
error analysis. The paper ends with some remarks in section |5] 



2 Preliminaries 



2.1 Finite difference approximations 

The spatial discretization is given by a set of points xj; j = 0, N and the discretization in time is 
represented by a sequence of times = < ... < < ... < T. For the sake of simplicity we will 
assume constant step sizes 6x and 6t in space and time respectively. The discrete solution at a point will 
be represented by n" w u{t^, xj). The schemes consist in computing approximate values n" of solution 
to O on [n6t, (n + l)6t[x [j6x, (j + l)6x[ forn G N and j G N thanks to the following relation: 



+ F{u]_ 



Ua , u 



u 



2< + < 1 



0, 



(9) 



where and F are, respectively, the discretizations of the nonlocal and nonlinear terms. Note that the 
Laplacian term is discretized using centred finite difference approximation. We begin by considering two 
discretizations l^^,!]^ for the operator X corresponding to formulae Q and respectively. In both 
cases, we use a basic quadrature rule on the mesh {[j6x, (j + l)Sx[)j^j^ to approximate each integral 
and we use a finite difference approximation of the derivative: 



5x-^/'Yl 
1=1 



■l+l 



(10) 



(11) 



1=1 ^ 

Let us remark that we begin the sums at / = 1 in order to avoid the singularity of l/|2|^/3 and l/|2;|^/3 
at z = 0. We will comment later on the truncation of the series, see Section HI Let us simply note that 
if (fj = for all j < then the series (ITOl ) is in fact a finite sum. Since the spatial mesh is given by 
{[j6x, {j + we will indeed assume that tpj = for all j < 0. 

Remark 2. Using fractional calculus, we could also consider, for any causal function if, the standard 
Griinwald-Letnikov formula for the fractional derivative X. Indeed, using the expression X can be 
approximated by the following two formulae 



r(2/3) 



l>0 



4/3 



r(2/3) 

6x^/'^ 



r(/ - 4/3) 



^r(/ + i)r( 

l>0 ^ ' ^ 



-4/3) 



and 



r(2/3) 



l>0 



-2/3 
I 



i^j-i+i - '^^j-i + ^j-i- 



(12) 



(13) 



where, for all a > and k £ N we denote by {^J binomial coefficient defined by 

a\ a{a — 1) ... (a — A; + 1) T{k — a) 

k) k\ " ^~ ' r(-a)r(A; + l) 
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and 



denotes the negative binomial given by 

_ p{p + 1) ■ ■ ■ {p + k - 1) 



(-1)' 



-p 



For more details about Griinwald-Letnikov derivatives, we refer the reader to the book U4\l . 

To analyze the stability of the discrete problem (|9]l using Fourier analysis, we investigate the follow- 
ing Unearized explicit scheme 



6t 

where w is a positive constant. 



6x 



+ r]Isx[u'" 



0, 



(14) 



Remark 3. In the case where we consider that v is a non-positive constant, dxU is discretized using a 
downstream finite difference approximation and so F is given by 



F{uU,u],u]^,) 



Therefore, taking into account the discretization ([TOl i. the numerical scheme is written as follows: 



,n+l 



1 2^^^ 
6x 6x'^ I 



V 6t e5t 



T] 6t 



+ 00 



(15) 



1=1 



and since 

+00 



1/3 



^j-l 



1=1 



1=2 

- u] -{2-2 



the numerical scheme ([TST i reads 

+1 _ 1^ n ,(i_!^_2l^-Jl^\u^^+f!!^ + l^ + (2-2-^/^)^\u'^ 



(5a;2 fe4/3 / 



(5a; (5x2 



(52;4/3 ; 



i-oo 

^[(z + i)-V3_2ri/3 + (^_i)-i 

i=2 



/3 



(16) 



Considering now the discretization ([TTI i. the numerical scheme (fT4l ) can be written as follows: 

e5t 



,n+l 



A r]6t 



V 6t e6t ^ V 6t eJt 



+ 00 



9 (5x^/3 



u 
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Recall that the Riemann zeta function, for Re(s) > 1 

oo 

as) = E 



n 



n=l 



Since 



+00 
1=1 



7/3 



u 



■j-i 



u 



2C(3)«"+i-C(5)«" 



)+oo 
1=2 



withcd: 



3.601, C(|) ~ 1.415, the numerical scheme reads 



u 



n+l 



e6t 



4 v^t 



6x^ 9^4/3 2^^3- 



+ 



v6t e6t A rjdt ,1^,4, , , „ 



4^^ 

9 (^a;4/3 ^ 



(17) 



Remark 4. If the Fowler equation ([T]) satisfied the maximum principle, a classical way to get sufficient 
conditions for the L"^ -stability of the scheme would be to ensure that li"^^ is a convex combination of 
{u^)j^^. Though one can easily check that all coefficients sum up to 1, we remark that {I + l)~^/3 _ 
2l~^/3 + — l)-i/3 > because theffinction x — > is convex and — | ^^fyg /~'^/3 ^ qj-qj- q// / > i 

Thus, uj~^^ is not a convex combination of To get conditions of numerical stability we have to 

rely on the Von Neumann method. 



2.2 The continuous amplification factor 

In this section, we are going to study the amplification factor of the following equation 

dtu(t,x) + vdxu{t,x) -edl^u(t,x) + r]I[u{t, -Mx) = 0. (18) 
Then, u(t, x) = e^^^-^"^ is a solution to ([T8] ) if and only if the following dispersion relation is satisfied 

(j + wk + ek^ - r?|A;|^/^^r(^) (l - iV3sign{k)) = 0, 

where A; G M and o" G C. Indeed, we have ut{t,x) = au{t, x),Ux{t, x) = iku{t,x),Uxx{'t,x) = 
—k^u{t, x) and 



r+00 

-k^u{t,x) / \i\-^/^e-'^^di, 



-k'^u{t, x) 



p+00 r+00 
Jo Jo 



\'^'lnl)+k\k\'/'^nl) 



u(t,x), 
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where we have used Fresnel integrals. 

Hence the mukipUcative factor which enables to get the solution at the time tn+i from the solution at the 
time tn is 

G,ont{k) = e-'"^(''\ (19) 
where (j){k) = ek^ - r/ir(|) |A;|^/3 + i (^r/^r(|) k\k\^/^ + u/c) . Therefore 



Ii-i f, y, y, f, }f,_ y 




;.;4.-; -.{jr;, -.fv-,. 



.A 



Figure 1: Behaviour of Re (cp) for -q, e fixed, fep = (^T (|) 2^'^/^ threshold frequency. 
Figure [T] shows that the modulus of the continuous amplification factor during one time step is con- 



trolled by e"*"^*, with := - minRe = - Re (t){k*) = ^ {^F (|))^ ^, where 



3 13/ e 



3/2 



Thereby, the exact continuous amplification is maximum for frequency fc*, and its modulus is bigger 
than 1 only for frequencies in the range (0,A;o]- The magnitude of this amplification during one time 
step will also be proportional to 6t. Obviously, this phenomenon affects only the low frequencies in 
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the range (0, ko] and strongly depends on the choice of parameters t] and e. This is why the standard 
definitions of stabihty are not adapted for this model because they do not take into account the possibility 
of amplification of certain frequencies. 

And since high frequencies are usually responsible of numerical instabilities, we are going to focus our 
attention on the high frequencies which are quickly dampened in Fowler's continuous model in order to 
exhibit numerical stability conditions. 

3 Stability analysis 

The purpose of this section is to study the numerical stability of schemes introduced in the previous 
section and to exhibit stability criteria. We recall that the numerical stability enables to ensure that the 
difference between the approximated solution and the exact solution remains bounded for all T > with 
6x, 6t given. To get numerical stability criteria, we consider the Von Neumann or Fourier method. In 
this approach, we assume that the discrete solution is written in as a single Fourier mode 

nj" = n2e*'=^% (20) 

where /c G Z is the wave number. Injecting (l20l ) in the numerical scheme ([141) . we get 

u''+' = g{5x,6t,k)ul, (21) 

where g is the discrete amplification factor. In what follows, for simplicity, we denote indifferently 

g{6x, 6t, k) = g{6x, 6t, 6), where 6 = k6x. 

Remark 5. Note that due to the aliasing phenomenon it is enough to study the discrete amplification 
factor for 9 G [0, vr]. 

Following the previous discussion concerning the notion of numerical stability (see Section |2]i, we 
introduce the following definition: 

Definition 1. We say that a numerical scheme which approximates the linearized Fowler equation prob- 
lem is stable if the high frequencies are strongly stable that is to say: 

3 < 6*0 < TT such that MO G (6*0, vr], \g{5x, 6t, 6)\ < 1, 

where g is the discrete amplification factor. 

Lemma 1. Let a, 6 G M and d G M+. Then we have 

V6I G [0, 27r],|a + 6e-*^| < d if and only if a + \h\ < d and a - \b\ > -d. 
Proof. We can easily check this property, see Figure |2] 

■ 

Proposition 1. The finite difference scheme (1141 ) is stable in the sense of Definition\J}if 6x and 6t satisfy 
the following conditions: 
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Figure 2: Dashed circle (resp. continuous circle) is centred at a (resp. 0) and of radius \b\ (resp. d). 



V 5t e5t , I/O, riSt 

^ + 2;^ + (2-2-'/=,£^<l, 



and if moreover, die space-step 6x is smaU enough in order that 



riSt e6t o.On. 



(22) 



v6t e5t 4 A ,4, A riSt 



(24) 



9(C(3)-l + C(3)J^<2^sm(-), (25) 

where Oq designates the stability threshold frequency. 
Proof. ForXj^.. 

For the numerical scheme ([T6] l. the amplification factor is given by: 

,^ ^ ^, V 5t e5t , 71 6t fv5t , M-\^'ri^t\ ,7) 

= ,___2-^(l-cos«)--^+ - + (2-2-'/')-^ 



00 



1/3 



e-'^\ (26) 
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where 6 = k5x. Since, for all G N 

N 



[(/ + 1)^1/3 _ 2/-1/3 + _ i)-l/3j = (JV + l)-^/3 _ ^-1/3 _ 2-1/3 ^ 1^ 

/=2 



then 



1/3 



1=2 



1 - 2 



-1/3 



> 0. 



Thus, from (|26] |. to have (^t, ^)| < 1 it is sufficient to have 



v6t 
6x 



n.O.eSt T]6t f V 5t , M-^^ ri 5t \ 



< 



r]5t 
5^ 



(1-2-1/3)^ 



where we assume that 



Next from Lemma [T] 



r]6t 



(1_ 2^1/3) < I 



fe4/3 

is satisfied if and only if we have 



v6t e5t 



l-2:-^-4-£lsin^(-)-(3-2-V3) 



6x 5x'^ 
A sufficient condition is then 



rj 6t 

fe4/3 



> 



1 - (1 - 2" 



-l/3^ 



rj5t 



5x^/3 J' 



(27) 



(28) 



(29) 



(30) 



Let < 6*0 < vr. Then, for all 6 G (^O) ^r], condition (l30l ) can be rewritten as 

Therefore, the numerical scheme (IT4] | with the discretization Xj^ is stable in the sense of Definition [T] if 
the space and time steps 6t, 6x satisfy the following conditions 



2 ■ 2 / ^0 N e 

(1 - 2-1/3) ™> (y);;. 



(31) 
(32) 



Note that from condition (l32l ). we can see that hypothesis (|29l ) is satisfied. 
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ForX? . Injecting 



g-2{5x,5t,e) 



in ([TT] ). the amplification factor g2 associated to this scheme is given by: 



v6t e6t 4 7 r]6t 



+ 



v6t e6t 4 14, 



_ 4 r]St ^,-7/3 -iei 



l>2 



v6t eSt 4^,7, r]6t 



.4 4 r/,^t 



^3^ fe4/3 



+ 



u5i _ 4 r/'^t \ _ 4 7?Jt ,-7/3 
6x 9 5x^/^J 9fe4/3Z^ 



7;(5t eJt 0,0. 4 



C(^)-c4)cosi 



r]5t 



+ 



7^ + 9(^(^ 



-i0 _ Y^^-7/3 

9 (5a;4/3 ^ 

/=2 



(33) 



Since 



+00 



J]r7/3 = c(^)-i«o.4i5, 



from (l33T l. 1(72(52;, 6t, 



vSt .e6t .2, 



< 1 if 
4 



where we assume that 



C(^)-C(^)cos, 



i]5t 

fe4/3 



fe4/3 j 



-i6» 



< 



1-^ C( 



7? 5t 



C( 



7, 



1 



rj 5t 

6x^/^ 



< 1. 



From Lemma [T] we have that (l34l ) is satisfied if and only if 
1 



e5t . n.e, 4 
^^^^^(2) + 9 



C(^)-l + C(^)(l-cos, 



U(5t e6t 0,6*, 4 A, 7, ^,4,, 



7? 5t 



< 1 



fe4/3 



> -1 + 



1 



fe4/3 ' 

(35) 

T] 6t 



(34) 



A sufficient condition is then 

4 



7, 



3' 
(5x 



-i + C( 

+ 2 



4 
3 

4 

fe2 + 9 



77 (^t e(5t 9,0, 



1 



'3' V (^x4/3 

Let < 6*0 < vr. Then, for all 6 G (^o, condition (l36l ) is rewritten as 



< 1. 



5 H' + 



(^x4/3 
7] St 
(5x4/3 



(36) 
(37) 

(38) 
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where C(i) - 1 + C(|) ~ 4.02. 

Note again that from condition (|37] ). we can see that hypothesis (1351) is satisfied. 



Notations. We will denote by CFL^^^ and CFL'^^^^ the following modified Courant-Friedrichs-Lewy 
conditions 

Some remarks. 

1. Condition (l22l ) (resp. (|23] )) can be seen as an extension of the classical CFL condition with in addition 
the anti-diffusive term jj^- This criterion is not more restrictive than the usual condition of stability 
without the nonlocal operator which corresponds to the linearized Burgers equation with viscous term. 
This condition is very restrictive on the space and time steps in particular because of the term |^ which 
stems from the explicit discretization of the Laplacian. In order to have less restrictive conditions, we 
can implicit some terms. For instance, if we decide to implicit the nonlocal and the Laplacian terms, 
condition (l22l) (resp. (1231 ) ) is reduced to 

v5t ^ 

6x 

We find again the well-known CFL condition. 

Figure [3] shows the behaviour of amplification factors for Zj^ and We can see, for Xj^, that the 
maximal value of 6t which ensures the numerical stability is 5tmax ~ 0.042 and that for this value we 
have CFLI^^^ 0.99. Figured displays the behaviour of the modulus of the amplification factor with 
discretization Xj^, as a function of 9. We can notice that the high frequencies are strongly amplified. This 
phenomenon illustrates the numerical instability because high frequencies should be quickly dampened. 

Figure |5] shows that the low frequencies are slowly amplified. This phenomenon is not due to the 
instability of numerical schemes but stems from the model. In Tables [T]|2] and [3] (see Section|4ll, we have 
studied the quotient t7=M— r , i = 1,2. We can see that globally the discrete schemes dampen more than 
the continuous problem when the stability conditions (1221) and (1231 ) are satisfied. 

2. Conditions (l24l) and (l25l) are unusual and deserve some explanations. The temi proportional to ^^^^3 
represents the amount of nonlocal anti-diffusion while the term proportional to |^ corresponds to the 
amount of classical diffusion. Both conditions simply mean that, for frequencies above the threshold 6*0, 
diffusion should control nonlocal anti-diffusion. 

We can see that conditions (l24l ) and (|25] ) cannot be satisfied for low frequencies. Indeed, for Oq close to 
0, these criteria impose to the space step to vanish, which is not possible. Let us note that this is coher- 
ent because the low frequencies are not "strongly stable", they are slowly amplified by the continuous 
problem. We can see in Figure |6]that condition (|24] | is not necessary. Indeed, if we choose the threshold 
Oq = tt/2 "large enough", condition (1241) reads 

6x < 0.25, (39) 
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2- I I I I 

Figure 3: Amplification factors for ij^ (blue line) and (dashed line). 




Figure 4: Amplification factor for Z^^ with CFL^^^^ ss 1.22. 

and we have plotted \gi \ in function of for 5x = 0.5 which does not satisfy the condition (|39l ) but we 
can still notice that the numerical scheme is stable. All numerical simulations that we performed confimi 
this statement. This leads us to think that condition (l24l ) (resp. (|25] )) is too pessimistic. In fact to estimate 
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Figure 5: Amplification factor Xj^ with t] = 8,v = l,e = 0.5 and 6x = 0.05, 6t = 0.001. For these 
coefficients, CFL]^^^ ^ 0.94. 



the magnitude of sums 



X;[(^ + i)"'/'-2ri/3 + (/- 1)^1/3 

1=2 

(resp. ^~^^^e~*'^), we just controlled the sum of the modulus 

oo 



-iW 



1=2 



(resp. ES^"^^^)- In this manner, we probably miss some cancellation effect of the e But we 
could not find any other way to estimate these polylogarithm series. 



3. Finally, in practice, the single condition (|22] | (resp. (l23T l) can be used to ensure the numerical stability 
of the scheme (fT4l ) withXj^ (resp. X|^). We saw in Figures |5]and|6]that the scheme with the discretization 
Xj^ is stable if condition (l22l) is satisfied. Figure [T] shows that the high frequencies are amplified, when 
condition (|23] ) is violated. This phenomenon is only due to numerical instability because the continuous 
problem quickly dampens the high frequencies. 
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Figure 6: Amplification factor gi for = u = 1, e = 0.1 and 6x = 0.5, 6t = 0.01. For these coefficients, 
we have CFLl, w 0.0584. 




Figure 7: Gain function for for r/ = u = e = 1, (5x = 0.1, (5i = 0.05. For these coefficients, we have 
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4 Error analysis 



4.1 Truncation error 

In this section, we analyze the truncation error. Finite difference scheme ([14] ) is consistent with the 
linearized partial differential equation if for any smooth function (p^t, x) the local error E^t^Sx satisfies 



Est,5x =P<i>-Pi 



5t,5x^ 



(40) 



as 5t, (5x — )■ with 



5t 



+ v 



^1 



' = (t>t + V (px - e 4>xx + ill[4>\, 



5x 



for i = 1,2. 



Remark 6. The practical implementation of the schemes requires to make some truncations. First, 
we consider a bounded domain [0, T] x [0, D] and to simplify, we also assume that 5t = T/Nst and 
6x = D/Nsxfor some integers Nsx and Nst- Another truncation concerns the integral operator for 
the nonlocal term I. We replace /q^°° with and in the finite difference approximations (1101 ) and 
(HD series Xy/^i '^^^ replaced with partial sums where A — Ag^ 5x. However, the truncation 

parameter A has to be chosen judiciously. A "short memory " principle has been investigated to choose 
this parameter. This principle is based on the fact that terms /~^/^ and in discretizations (1101 ) and 
(1111 ) decrease with I therefore, we have to take into account the behaviour of (p{x) only in the recent 
past, i.e. in the interval [x — L,x], where L > is called the "memory length". Finally, the use of the 
short-memory principle leads to the simple replacement of^^^ by 'Ylt=i' where A^^ = [^] 
Note that the truncation parameter A also strongly depends on the discretization of the nonlocal term 
X because /^^/^ decreases more quickly than For the sake of simplicity, we will denote by A the 

truncation parameter for the discretizations dlOb and dill) . 

Proposition 2 (Local eiTor). Tlie local error of the numerical scheme (1141) satisfies: 



• Porll: 

\Est,sx,A\ ^ ^('^*) + 0('^a;2/3) + 0(^-i/3) + Q (^^-4/3 ^ ^^2 ^^3^ + Q {A fe^) . (42) 



Proof. From Taylor series, we have 



LTl + l 



(t>t {tn, Xi) - = 0{6t), (43) 

or 
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<t^x {tn, Xi) - /'-^ = 0[bx), (44) 

OX 

(pxx [tn ,Xi) ^ =0[dx^). (45) 



Let us now study the truncation error for the nonlocal term X. 
For Xj^ : We rewrite dD as follows 



.•1 J(i-l/2)5x Jo 



'(j-l/2)fc 70 

r-A f+OC 

(46) 



/"A r+oo 



and the discretization (ITOl ) becomes 



'(i-l/2)5x 



with := j(52; and = -^-2+1 — — i;^?^ Using (l46l ). we then get the following relation: 

A^:, „(j+l/2)5x 



+00 



= Ti-T2+T3-T^. 

Let us study the term Ti. Since 

^Jx^''^sx-r'^'^xx{t'',x,-o = (Q,'^'-r'^')^xxir,x,-o 



then 



+ ^Jx^'i^'sx-4>xx{t'',X,-0) 



Asx p{j+l/2)Sx 

E / i^Jx^' - r'/' ) </'xx(i", Xi - 

P^J{j~l/2)5x ^ ' 
Mx /.(j+l/2)<5a; 

+E / Cx'\nx-'^ut''.^i-i))di. 
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By the mean value theorem appUed to z — )• \z\ we have for all ^ G [{j — \)5x] (j + \)5x\ 



-1/3 



Sx 



e-l/3| < 



sup 

z£[(j~\)5x;{j + \)5x] 



< l\{j-\)H^'^'\Csx-Cl 

< ^\ij-^)5x\~*/Hx. 



Thus, integrating over [{j — ^)Sx; (j + ^)Sx] we get 



(47) 



because X]j>i (j-i/2)4/a ^ ^'^^ ^ a positive constant which depends on | |i;^>xx| |l°o((o,t)xR)- 
Moreover, by classical midpoint quadrature rule 

l'{j+l/2)6x 

/ (pxx 
J(j-l/2)5x 

with r]j G [{j - l/2)6x; (j + l/2)6x] then 



{t'^,Xi - ^) = 6x(j)xxit"',Xi - (,sx) + -^4'4x{t"',Xi - rjj), 



Ti,2 = Y.(j5x) 

= ^(jSx) 
From Taylor series, we have 



-1/3 



-1/3 



^^'^Sx 



u+imsx 



ij-l/2)5x 
fe2 



^(f'ixi't^ , Xi — r]j) 



ct>xx{t^,x^ - isx) = - ^'A4x(i",2;i - isx) + 0{5x^), 



thus we obtain 
Ti,2 ^ 
We finally get 



As 

1 



^U^x) 



i=2 



which imphes that 



(48) 



(49) 



(50) 
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because Asx = 

We next control the term T2 by 

IT2I < ii</'..iil-((o,t)xr) / ' ier'/'de = cfe2/3. (sd 

We estimate as follows: 

i-A+Sx 

\T'i\ < ||<^xx||l°°((o,t)xR) / di, 

J A 

< \\4>xx\\l^{{o,t)xM.)A~^^^^x. (52) 

Finally, using an integration by parts, the term T4 is written as 

r+00 

J A 

1 r+oo 

hence, we obtain 

IT4I < CA-^/^ (53) 

where C is a positive constant which depends on ||i?!*a;||L°°((o,T)x]R)- 

Hence, using relations (03]), (05]), (07]), dSO]), dSB, ^ and (ES), we obtain 

\El,st,A\ = \P^itn:X^)-Plsx^\ < 0(fe2/3) + 0(5t) + 0(A-V3) 

which completes the proof for Xj^. 
For : As previously, we rewrite (jS) as 



J{j-l/2)5x Jo 
pA r+oo 



with <I>(t", = I Xi-^) - (/>(t", Xj) + (/>a;(t", Xi) and the approximated integral (HB be- 

comes 

fa p(j+l/2)Sx 



^(j-l/2)5x 

with = I (</.f_ . + ^^ii^i) and = j5x. 



19 



Let us now estimate the error on the nonlocal term. 

„(j+l/2)<5a: 



,4 



+00 



/■+00 

J A 

= Ti-T2 + n-T^. 
Let us study the term Ti. As previously for Xj^., we rewrite Ti as 

^"5^ f(i+l/2)<5x 
j.^l ^(j-l/2)5x 

By the mean value theorem applied to z — )• we have for all ^ G [(j — l/2)fe; (j + l/2)(5a 



6 |(j-l/2)5x|io/3' 



6 (j - 1/2)10/3 

Next, by Taylor-Lagrange formula, we have 

4 ^ r{j+l/2)5x 
'(i-l/2)5x 



9 Jij~l/2)5x 



^■=1 V(j-l/2)5x 

^fc -I /■(7+l/2)5x 

U - l/2)io/'^ J(j~l/2)Sx 

■ 2 

because X^jLi (j^i/2)"V3 < co. C denotes a positive constant which depends on | |l°°({o,t)xM) ^rid 
may vary from line to line. 

Moreover, using again midpoint quadrature rule, we have 



/.(j+l/2)<5x- r 3 

J{j-l/2)5x ^4 
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with rjj G [(j — ^)(5x, {j + ^)6x]. Hence, 



Tl,2 = ^ Csx 



7/3 



6x^ 



But using again Taylor expansion, we get 



2dx D 



and so 
Thus, 

^1,2 = 



" + , j - ./-(t", - jfe) + <P{e, Xi) - x,)j6x 



Sx' 



Csx"-^cPs,{t'',Xi)+jO{6x^) 





^5x 



6x^ 



6x^ 



7/3 



6x^ 



4 6x^ 



^Sx—hxit'',Xi)+jO{6x'') - -^(t)xx{t^,x^ - r]j) 



We finally get 



^5x 



= 0{6x^^''^) + 0{A^Sx^) +0{A6x^) , 

where C is a positive constant which depends on | |i;^3x| |l°°((o,t)xR) II</'xx||l°°((o,t)> 
Let us now study T2. Using Taylor-Lagrange formula, we have 



\T2\ < C\\(j)xx\\L°°{{0,T)> 



— a2 
2 



ie|V3 



where C is a positive constant. 
Let us next consider T^: 
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Then 



m < c [ " m-'^' + m-^/')dt 

J A 

< C (^A~'/^ + A-^/^) 6x, 

< O (^5x A-^/'^y 

with C apositive constant which depends on ||</'||l°°((o,t)xIR) and ||(/>x||l°°((o,t)xIR)- 
And since 

A r+oo A r+oo 

< CA-^/^, 

where C is apositive constant which depends on | If/*! |l°o((o,t)xM) and | |(/>x| |l°o((o,t)xR)' we finally get 

+ 0{A^5x'^) +0{A5x^) . 
The proof of this proposition is now completed. ■ 

Remark 7. From previous Proposition, we can see that the numerical scheme (114b with Xj^ ( resp. j 
is consistent if 5x « A^'^^'^ (resp. 5x « A^"^!^). 

4.2 Convergence experiments. 

In this section, we investigate the convergence using numerical simulations. Despite much effort we are 
unable to prove theoretically the convergence of the numerical solution towards the exact continuous 
solution. Indeed, the Lax procedure "stability + consistence = convergence" cannot be applied here due 
to the instability of low frequencies. 

In what follows, Z^-norm is used to measure the accuracy of approximated solutions. Thus, we analyze 
the following error 

1 ^ 

Ex = j^Y.^\u](T)-u]{T)\\ (54) 

n=0 

where u^, -u^ are, respectively, computed for space steps (5x/2 and bxjA:, until a final time T. 
Figure |9] shows the numerical convergence rates obtained with the initial data displayed in Figure [8] 
These rates were obtained using bx = 10~^, 10~^, 10""^, 10~^. We plot the logarithm of the eiTor Ei 
versus the logarithm of 5x. 

4.3 Phase error 

Numerical schemes produce, according to cases, results ahead of or delayed w.r.t exact solutions. In this 
section, we are interested in the error made on the velocity introduced by the discretization. Let us first 
note that, in addition to the anti-diffusive effect, the nonlocal term is also responsible of the motion of 
the initial data. Indeed, we saw in Section |2] that the continuous amplification factor has an imaginary 
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Figure 8: Initial data used for numerical experiments. 




lea . I't ^1 



Figure 9: Convergence in space for Xj^ (red) and I^^ (blue). Dotted line has slope 2/3. 

-.j/'n v^r(-) k\k\^^-^+ vk] St 

part e V 2 ^ 3 ^ I I y Therefore, the advection temi is not the unique factor of displacement. 
It is the error on the argument —6t{v k + ■^r(|)?7 A;|fc|^/'^) that causes the phase error. To evaluate this 
error, we rewrite the discrete amplification factor gj introduced in Section |3] as 

9j = \9j\e ^ 

for j = 1,2, where 6^. is the argument of the discrete ampUfication factor gj. The phase lag during one 
time step is then given by 

Ej = {vk + ^T{^)r^k\k\'/^)5t-0,^. 
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Thus, if Ej is positive, the numerical wave goes slower than the physical wave and it goes faster if Ej is 
negative. We have computed, in Tables [T] [2] and [3l the phase delay after one oscillation 

Aj = l , for j = 1,2, 

(^yk + ^T{l)r]k\k\y^)5t 

for different values of Cr = ^,Df = ^ and Fo = j^. 

We remark that scheme with the discretization involves a delay larger than the model with Xj^. 

4.4 Discrete amplification vs. continuous amplification factors 

Let us define 

r< |gil n - l^2| 



I cont I I ^cont \ 

for different values of Cr = Dj = ^ and Fo = j^. Results are reported in Tables [I] |2] and 
111 We have CFL]^^^ = Cr + Df + XiF^, where Ai = 2 - 2-1/3, Aa = | (C(|) - l) • When the 
condition CFL\^^^ > 1 is violated because Cr or Df are close to one, we note that high frequencies 
9 > 7r/2 are more amplified by the discrete schemes than by the exact continuous problem. Whereas 
when the condition CFU^^^ > 1 is violated because Fq is close to one, high frequencies 9 > tt/2 
we, less amplified by the discrete schemes than by the continuous problem, as can be checked in Table 
121 This is one unexpected benefit of the discretization scheme in the unfavourable case where nonlocal 
anti-diffusion is predominant. If we take a closer look at Table [3l we notice that \g2 \ may be greater than 
one even if CFL"^^^ < 1. This is due to the fact that rj being big, the stability threshold frequency is 
close to vr, the aliasing limit frequency. 



Cr 






9 


Ai 


A2 


Gi 


G2 


0.2 


0.5206 


0.5156 


7r/6 


0.0082 


0.0333 


0.9584 


0.9788 








7r/4 


0.0024 


0.0573 


0.9102 


0.9550 








7r/2 


-0.0715 


0.1610 


0.6824 


0.8394 








37r/4 


-0.2684 


0.3911 


0.3788 


0.6626 


0.5 


0.8206 


0.8156 


7r/6 


-0.0128 


0.0104 


0.9541 


0.9736 








7r/4 


-0.0433 


0.0091 


0.9048 


0.9452 








7r/2 


-0.02476 


-0.0315 


0.7439 


0.8152 








37r/4 


-0.4733 


-0.1628 


0.8284 


0.6391 


0.9 


1.2206 


1.2156 


7r/6 


-0.0103 


0.0107 


0.9870 


1.0047 








7r/4 


-0.0306 


0.0128 


0.9869 


1.0182 








7r/2 


-0.0781 


0.0172 


1.1776 


1.1522 








37r/4 


-0.0210 


0.0371 


1.8187 


1.5053 



Table 1: Dampening and phase eiTor for -D/ = 0.2 and = 0.1 
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e 


Ai 


A2 


Gi 


G2 


0.2 


0.4206 


0.4156 




0.0213 


0.0481 


0.9654 


0.9859 








1 A 

'k/A 


0.0303 


0.0868 


0.9250 


0.9707 








7r/2 


0.0525 


0.2562 


0.7340 


0.9057 








O 1 A 

37r/4 


0.1721 


0.5567 


0.4834 


0.8487 


0.4 


0.6206 


0.6156 


7r/6 


-0.0066 


0.0216 


0.9649 


0.9860 








7r/4 


-0.0358 


0.0276 


0.9216 


0.9701 








7r/2 


-0.3470 


0.0154 


0.6750 


0.8841 








37r/4 


-2.0019 


0.0277 


0.4153 


0.7059 


0.8 


1.0206 


1.0156 


7r/6 


-0.0673 


-0.0361 


0.9614 


0.9837 








7r/4 


-0.1989 


-0.1169 


0.9022 


0.9568 








7r/2 


-3.1203 


-1.4458 


0.5305 


0.6566 








37r/4 


-3.8370 


-3.6360 


5.5254 


3.5099 



Table 2: Dampening and phase enor for Cr = 0.1 and Fq = 0.1 



Fo 






e 


Ai 


A2 


Gi 


G2 


Ifil 


\92\ 


Gcont 


0.2 


0.5413 


0.5312 


7r/6 


0.0307 


0.0784 


0.9455 


0.9852 


0.9741 


1.0150 


1.0302 








7r/4 


0.0363 


0.1360 


0.8852 


0.9707 


0.9180 


1.007 


1.0371 








7r/2 


-0.0079 


0.3495 


0.6236 


0.9052 


0.6240 


0.9057 


1.0005 








37r/4 


-0.1876 


0.6374 


0.3216 


0.8182 


0.2822 


0.7180 


0.8776 








vr 


-1.2548 


1.0000 


0.0574 


0.6017 


0.0574 


0.6017 


0.6950 


0.5 


0.9031 


0.8780 


7r/6 


0.0591 


0.1552 


0.8992 


0.9875 


1.0092 


1.1084 


1.1224 








7r/4 


0.0650 


0.2538 


0.8043 


0.9758 


0.9664 


1.1725 


1.2016 








7r/2 


-0.0103 


0.5264 


0.5235 


0.8684 


0.7589 


1.2590 


1.4498 








37r/4 


-0.0906 


0.7619 


0.4215 


0.6524 


0.6992 


1.0824 


1.6590 








vr 


-0.0430 


1.0000 


0.4202 


0.5110 


0.7435 


0.9042 


1.7694 


0.9 


1.3857 


1.3404 


7r/6 


0.1064 


0.2438 


0.8602 


0.9941 


1.0822 


1.2512 


1.2583 








7r/4 


0.1361 


0.3752 


0.7523 


0.9762 


1.0999 


1.4273 


1.4621 








7r/2 


0.2002 


0.6556 


0.5123 


0.7449 


1.2178 


1.7707 


2.3771 








37r/4 


0.2981 


0.8366 


0.3873 


0.4083 


1.5020 


1.5836 


3.8781 








vr 


0.3924 


1.0000 


0.2696 


0.2126 


1.6583 


1.3075 


6.1519 



Table 3: Dampening and phase error for Cr = 0.1 and Df = 0.2. 
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5 Concluding remarks 



We have presented in this work a first investigation of finite differences schemes approximating the 
Fowler equation. We saw that the anti-diffusive behaviour of the nonlocal term does not enable to con- 
sider the classical notion of stability. Nevertheless, considering only the behaviour of the high frequen- 
cies (which should be quickly dampened), we exhibit numerical stability criteria which can be used to 
make simulations. Numerical computations have shown that numerical schemes dampened more than 
the continuous problem. Finally, consistency property has been proved and convergence of schemes has 
been investigated. 
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